#delimit;
version 13;
set more off;
program drop _all;
  quietly log;
  local logon = r(status);
  if "`logon'" == "on" {; log close; };
log using bos2015prq-robust.log, text replace;


/* *************************************************************      */
/*  File Name:  bos2016prq-robust.do                                  */
/*  Authors:    Frederick J. Boehmke and Emily Schilling              */
/*  Date:       May 03, 2018                                          */
/*  Purpose:    Replication of pivotal politics and initiative use    */
/*              models for Boehmke, Osborn, and Schilling, SPPQ,      */
/*		2015. This runs the robustness analysis reported on   */
/*              at the end of the section "Pivotal Politics and       */
/*              Initiative Use."                                      */
/*  Input:      bos2015prq.dta                                        */
/*              state legislator scores may 2013.dta                  */
/*  Output:     bos2015prq-robust.log,                                */
/*              bos2015prq-robust.gph                                 */
/* *************************************************************      */


	/************************************************/
	/* This program combine the estimates from the 	*/
	/* multiply imputed data sets.			*/
	/************************************************/


program define micalc, eclass;

  syntax, PARameters(name) COVariance(name) imputations(integer);

	local b `parameters';
	local C `covariance';
	local m = `imputations';
  
	local names: colnames(`b'1);

	forvalues i=1/`m' {;

	  if `i'==1 {;
		matrix b = b1;
		matrix C = C1;
		};

	  else {;
  		matrix b = b + b`i';
		matrix C = C + C`i';
		};
		
	  };

	matrix b = b/`m';
	matrix C = C/`m';
	
	forvalues i=1/`m' {;

	  matrix bb=b`i'-b;

	  if `i'==1 {;
		matrix Vb=(bb)'*bb;
		};

	  else matrix Vb = Vb + (bb)'*bb;

	  };

	matrix Vb=Vb/(`m'-1);				/* Covariance adjustment */

	matrix V=C+(1+1/`m')*Vb;				/* Final covariance   */

	ereturn post b V;
	ereturn display, plus neq(3);
	

end;


use bos2015prq, clear;


	/********************************/
	/* This is the reference model. */						
	/********************************/

	
nbreg init_use2_iri filibint_sm vetoint_sm statinit consinit init_dist init_sigs init_circ init_circ_unl 
	unif_any totpop rpcpinc citi6010 if mod(year,2)==0, robust;

	
	/******************************************************************/
	/* Read in the Shor and McCarty data on ideal points. 		  */						
	/******************************************************************/

	
tempfile shormccarty;
	
use "state legislator scores may 2013.dta", clear;

		/* Combine the two chambers into one. */

  reshape long senate house, i(st_id st) j(year);

  generat state=st;

  replace np_score = . if senate == . & house == .;
  
  save `shormccarty', replace;
  
		/* Calculate the median legislator in each states. */
  
	collapse (p50) median_sm=np_score, by(state year);
	
		/* Merge back into the data with all legislators ideal points. */
  
	merge 1:m state year using `shormccarty', assert(match);
	
	drop _merge;
	
	drop if missing(np_score);
	
  save `shormccarty', replace;
	
	
	/******************************************************************/
	/* Assume median voter is a linear function of median legislator. */						
	/******************************************************************/

clear;

		/* Set number of observations based on the number of parameter configurations to run. The parameters */
		/* correspond to the representation function legislator = a + b*district - e*{Dem} + e{Rep} + u. */
	
set obs `=9*10*4*5';
		
generat a = .;
generat b = .;
generat e = .;
generat sd = .;
	
		/* Create variables in which to store results. */
	
generat sig_filint_10 = 0;
generat sig_vetint_10 = 0;

generat sig_filint_05 = 0;
generat sig_vetint_05 = 0;

generat sig_filint_01 = 0;
generat sig_vetint_01 = 0;

generat b_filint = .;
generat b_vetint = .;

generat se_filint = .;
generat se_vetint = .;

		/* Save the data set in which to start accumulating the results. */

save bos2015prq-robust, replace emptyok;

		/* i will index the iteration number. */

local i=1;

nois _dots 0, title(Loops running) reps(`=_N');

		/* Loop over all the parameter combinations. */

foreach a of numlist -0.2(.05)0.2 {;

  foreach b of numlist 0.4(0.1)1.3 {;
  
	foreach e of numlist 0(0.1)0.3 {;
	
	  foreach sd of numlist 0(0.05)0.2 {;

		forvalues j = 1/10 {;
		  
		  quietly {;
	    
			use `shormccarty', clear;
			
				/* Generate district median voters' ideal points based on the inverted representation function. */
				/* This assumes a linear relationship between the representation and the median voter. Here we */
				/* start with the representative's observed ideal point and back out the implied median voter ideal point. */ 
				/* Here is the representation function: np_score = a + b*district - e*{Dem} + e{Rep} + u. */
				/* Ideal points are then adjusted relative to the overall legislative median: */
				/* np_score - median_sm = a + b*(district-median_sm) - e*{Dem} + e{Rep} + u. */
			
			generat district = (1/`b')*((np_score - median_sm) - `a' - `e' - `sd'*invnorm(uniform())) + median_sm if party == "R";
			replace district = (1/`b')*((np_score - median_sm) - `a' + `e' - `sd'*invnorm(uniform())) + median_sm if party == "D";
			replace district = (1/`b')*((np_score - median_sm) - `a' +  0  - `sd'*invnorm(uniform())) + median_sm if party == "I";

				/* Then calculate the implied median voter in the state. */
			
			collapse (p50) med_voter=district, by(state year);
			
				/* Merge in the analysis data, create the implied intervals, and estimate the model. */
			
			merge 1:1 state year using bos2015prq, keep(match) generat(_merge_sim);

			generat filibint_sm_sim = abs(filib_sm - med_voter);
			generat vetoint_sm_sim = abs(veto_sm - med_voter);

			generat inside = .;
			replace inside = 1 if veto_sm  <= med_voter & med_voter <= filib_sm;
			replace inside = 1 if filib_sm <= med_voter & med_voter <= veto_sm;

			nbreg init_use2_iri filibint_sm_sim vetoint_sm_sim statinit consinit init_dist init_sigs init_circ init_circ_unl 
			  unif_any totpop rpcpinc citi6010 if mod(year,2)==0, robust iterate(100);
			  
			  matrix b`j' = e(b);
			  matrix C`j' = e(V);
			  
			use bos2015prq-robust, clear;
			
				/* Increase the tally of significant results if the model */
				/* converged within the 100 iterations maximum.           */
			
			if e(ic) < 100 {;
			  
				replace sig_filint_10 = sig_filint_10 + 1 if _b[filibint_sm_sim]/_se[filibint_sm_sim] > 1.65 in `i';
				replace sig_vetint_10 = sig_vetint_10 + 1 if _b[vetoint_sm_sim]/_se[vetoint_sm_sim] > 1.65 in `i';
				  
				replace sig_filint_05 = sig_filint_05 + 1 if _b[filibint_sm_sim]/_se[filibint_sm_sim] > 1.96 in `i';
				replace sig_vetint_05 = sig_vetint_05 + 1 if _b[vetoint_sm_sim]/_se[vetoint_sm_sim] > 1.96 in `i';
				  
				replace sig_filint_01 = sig_filint_01 + 1 if _b[filibint_sm_sim]/_se[filibint_sm_sim] > 2.58 in `i';
				replace sig_vetint_01 = sig_vetint_01 + 1 if _b[vetoint_sm_sim]/_se[vetoint_sm_sim] > 2.58 in `i';
				
				};
			
			save bos2015prq-robust, replace;

			};

		  };
		  
		  use bos2015prq-robust, clear;
			
			/* Calculate the average model results across the 10 draws and save the coefficients of interest. */
			
		  quietly {;
			  
			micalc, parameters(b) covariance(C) imputations(10);

			replace b_filint = _b[filibint_sm_sim] in `i';
			replace b_vetint = _b[vetoint_sm_sim] in `i';

			replace se_filint = _se[filibint_sm_sim] in `i';
			replace se_vetint = _se[vetoint_sm_sim] in `i';

			replace a  = `a'  in `i';
			replace b  = `b'  in `i';
			replace e  = `e'  in `i';
			replace sd = `sd' in `i';
			
			save bos2015prq-robust, replace;

			};
		
		  nois _dots `i' 0;
		  
		  local i = `i'+1;
		  
		};  
		  
		
	  };
	 	
	};
	
  };

		/* Create a table of the proportion of significant results across the parameter combinations. */
  
bysort sd: table a b e, c(mean sig_filint_10);
bysort sd: table a b e, c(mean sig_vetint_10);
  
bysort sd: table a b e, c(mean sig_filint_05);
bysort sd: table a b e, c(mean sig_vetint_05);
  
bysort sd: table a b e, c(mean sig_filint_01);
bysort sd: table a b e, c(mean sig_vetint_01);
  
label variable a 	"Intercept";
label variable b 	"Slope";
label variable sd 	"Standard Deviation";
  
		/* Create string variables to facilitate the graph below. */
  
generat e_str = "0" if e==0;
replace e_str = "0.1" if abs(e-0.1)<0.01;
replace e_str = "0.2" if abs(e-0.2)<0.01;
replace e_str = "0.3" if abs(e-0.3)<0.01;
  
generat sd_str = "0" if sd==0;
replace sd_str = "0.05" if abs(sd-0.05)<0.01;
replace sd_str = "0.10" if abs(sd-0.10)<0.01;
replace sd_str = "0.15" if abs(sd-0.15)<0.01;
replace sd_str = "0.20" if abs(sd-0.20)<0.01;
  
generat zero = 0;
generat one = 1;

		/* Calculate measures for positive/none/negative significance and magnitude (-3 to 3). */

generat sig_filint = sign(b_filint) if abs(b_filint/se_filint) > 1.65;
generat sig_vetint = sign(b_vetint) if abs(b_vetint/se_vetint) > 1.65;
  
replace sig_filint = 2*sign(b_filint) if abs(b_filint/se_filint) > 1.96;
replace sig_vetint = 2*sign(b_vetint) if abs(b_vetint/se_vetint) > 1.96;
  
replace sig_filint = 3*sign(b_filint) if abs(b_filint/se_filint) > 2.58;
replace sig_vetint = 3*sign(b_vetint) if abs(b_vetint/se_vetint) > 2.58;
  
keep a-sig_vetint;
  
compress;

save bos2015prq-robust, replace;

  
	/*******************************************************/
	/* Now graph the results of the robustness simulation. */
	/* These are basically heat maps of the proprtion of   */
	/* cases that are positive and significant.            */
	/*******************************************************/

	
foreach sig in 01 05 10 {;
  
	twoway scatter a b if sig_filint_`sig'==10, by(e_str sd_str, note("") legend(off) title(Filibuster Interval)) scheme(s1mono)
		msymbol(S) msize(*3.25) mcolor(gs0)
	  ||	scatter a b if sig_filint_`sig'==9,
		msymbol(S) msize(*3.25) mcolor(gs2)
	  ||	scatter a b if sig_filint_`sig'==8,
		msymbol(S) msize(*3.25) mcolor(gs4)
	  ||	scatter a b if sig_filint_`sig'==7,
		msymbol(S) msize(*3.25) mcolor(gs6)
	  ||	scatter a b if sig_filint_`sig'==6,
		msymbol(S) msize(*3.25) mcolor(gs8)
	  ||	scatter a b if sig_filint_`sig'==5,
		msymbol(S) msize(*3.25) mcolor(gs10)
	  ||	scatter a b if sig_filint_`sig'==4,
		msymbol(S) msize(*3.25) mcolor(gs12)
	  ||	scatter a b if sig_filint_`sig'==3,
		msymbol(S) msize(*3.25) mcolor(gs13)
	  ||	scatter a b if sig_filint_`sig'==2,
		msymbol(S) msize(*3.25) mcolor(gs14)
	  ||	scatter a b if sig_filint_`sig'==1,
		msymbol(S) msize(*3.25) mcolor(gs15)
	  || 	rline one one a, horiz
		lcolor(red)
	  || 	rline zero zero b,
		lcolor(red)
		ylabel(#8, labsize(small)) xlabel(#10, labsize(small)) 
		xsize(6.5) ysize(6.5)
		aspectratio(0.9)
		name(filint, replace);
	  
	twoway scatter a b if sig_vetint_`sig'==10, by(e_str sd_str, note("") legend(off) title(Veto Interval)) scheme(s1mono)
		msymbol(S) msize(*3.25) mcolor(gs0)
	  ||	scatter a b if sig_vetint_`sig'==9,
		msymbol(S) msize(*3.25) mcolor(gs2)
	  ||	scatter a b if sig_vetint_`sig'==8,
		msymbol(S) msize(*3.25) mcolor(gs4)
	  ||	scatter a b if sig_vetint_`sig'==7,
		msymbol(S) msize(*3.25) mcolor(gs6)
	  ||	scatter a b if sig_vetint_`sig'==6,
		msymbol(S) msize(*3.25) mcolor(gs8)
	  ||	scatter a b if sig_vetint_`sig'==5,
		msymbol(S) msize(*3.25) mcolor(gs10)
	  ||	scatter a b if sig_vetint_`sig'==4,
		msymbol(S) msize(*3.25) mcolor(gs12)
	  ||	scatter a b if sig_vetint_`sig'==3,
		msymbol(S) msize(*3.25) mcolor(gs13)
	  ||	scatter a b if sig_vetint_`sig'==2,
		msymbol(S) msize(*3.25) mcolor(gs14)
	  ||	scatter a b if sig_vetint_`sig'==1,
		msymbol(S) msize(*3.25) mcolor(gs15)
	  || 	rline one one a, horiz
		lcolor(red)
	  || 	rline zero zero b,
		lcolor(red)
		ylabel(#8, labsize(small)) xlabel(#10, labsize(small)) 
		xsize(6.5) ysize(6.5)
		aspectratio(0.9)
		name(vetint, replace);
	  
	graph combine filint vetint, scheme(s1mono) 
		rows(1) imargin(tiny)
		xsize(6.5) ysize(3.25)
		l1title(Size of the partisan shift increases down rows) 
		b1title(Variance of the random error increases across columns)
		saving(bos2015prq-robust-sig-pct-alpha-`sig', replace);

	};
		
  
	/*********************************************************/
	/* Graph the significance averaged across the 10 errors. */
	/* These are basically heat maps of the level of         */
	/* significance (.01, .05, .1) for positive coefficients.*/						
	/*********************************************************/


bysort sd: table a b e, c(mean sig_filint);
bysort sd: table a b e, c(mean sig_vetint);
  
twoway scatter a b if sig_filint==3, by(e_str sd_str, note("") legend(off) title(Filibuster Interval)) scheme(s1mono)
	msymbol(S) msize(*3.25) mcolor(gs0)
  ||	scatter a b if sig_filint==2,
	msymbol(S) msize(*3.25) mcolor(gs8)
  ||	scatter a b if sig_filint==1,
	msymbol(S) msize(*3.25) mcolor(gs13)
  || 	rline one one a, horiz
	lcolor(red)
  || 	rline zero zero b,
	lcolor(red)
	ylabel(#8, labsize(small)) xlabel(#10, labsize(small)) 
	xsize(6.5) ysize(6.5)
	aspectratio(0.9)
	name(filint, replace);
  
twoway scatter a b if sig_vetint==3, by(e_str sd_str, note("") legend(off) title(Veto Interval)) scheme(s1mono)
	msymbol(S) msize(*3.25) mcolor(gs0)
  ||	scatter a b if sig_vetint==2,
	msymbol(S) msize(*3.25) mcolor(gs8)
  ||	scatter a b if sig_vetint==1,
	msymbol(S) msize(*3.25) mcolor(gs13)
  || 	rline one one a, horiz
	lcolor(red)
  || 	rline zero zero b,
	lcolor(red)
	ylabel(#8, labsize(small)) xlabel(#10, labsize(small)) 
	xsize(6.5) ysize(6.5)
	aspectratio(0.9)
	name(vetint, replace);
	
graph combine filint vetint, scheme(s1mono) 
	rows(1) imargin(tiny)
	xsize(6.5) ysize(3.25)
	l1title(Size of the partisan shift increases down rows) 
	b1title(Variance of the random error increases across columns)
	saving(bos2015prq-robust-sig-avg, replace);
  
  
log close;
clear;
exit, STATA;
